Multi Anova

libraries

library(foreign)
library(ggplot2)
library(plyr)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(car)

Loading required package: carData
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘car’

The following object is masked from ‘package:dplyr’:

    recode

library(haven)
library(sjstats)

Registered S3 methods overwritten by 'lme4':
  method                          from
  cooks.distance.influence.merMod car 
  influence.merMod                car 
  dfbeta.influence.merMod         car 
  dfbetas.influence.merMod        car 

library(rstatix)

Attaching package: ‘rstatix’

The following objects are masked from ‘package:plyr’:

    desc, mutate

The following object is masked from ‘package:stats’:

    filter

library(psych)

Attaching package: ‘psych’

The following object is masked from ‘package:sjstats’:

    phi

The following object is masked from ‘package:car’:

    logit

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

data <- read.dta("/Users/aydin/Desktop/stat_pro/Gruppe A/Aufgabe5/WVS_Cross-National_Wave_7_stata_v1_6.dta")

#View(data)

plyr::count(data, 'C_COW_ALPHA')

plyr::count(data, 'Q241')

plyr::count(data, 'Q242')

plyr::count(data, 'Q243')

plyr::count(data, 'Q244')

plyr::count(data, 'Q245')

plyr::count(data, 'Q246')

plyr::count(data, 'Q247')

plyr::count(data, 'Q248')

plyr::count(data, 'Q249')

Since numeric variables are required, string variables should be assigned numeric values

data['dem1'] <- mapvalues(data$Q241, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem2'] <- mapvalues(data$Q242, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem3'] <- mapvalues(data$Q243, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem4'] <- mapvalues(data$Q244, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem5'] <- mapvalues(data$Q245, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem6'] <- mapvalues(data$Q246, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem7'] <- mapvalues(data$Q247, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem8'] <- mapvalues(data$Q248, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))
data['dem9'] <- mapvalues(data$Q249, from = c('It is against democracy (spontaneous)','Not an essential characteristic of democracy','An essential characteristic of democracy'), to = c(0,1,10))

Create a new column as a dependent variable

data['dem_num']<- (as.numeric(data$dem3) - as.numeric(data$dem5) + as.numeric(data$dem6) - as.numeric(data$dem8) + as.numeric(data$dem9)) / 5

Histogram for dependent variable

hist(data$dem_num, main = "Histogram", ylab = "Count", xlab = "demoracy perception", col = "lightblue")

factors reducing affecting the level of education

data['edu'] <- mapvalues(data$Q275, from = c('Early childhood education (ISCED 0) / no education',
                                    'Primary education (ISCED 1)',
                                    "Lower secondary education (ISCED 2)",
                                    "Upper secondary education (ISCED 3)",
                                    "Post-secondary non-tertiary education (ISCED 4)",
                                    "Short-cycle tertiary education (ISCED 5)",
                                    "Bachelor or equivalent (ISCED 6)",
                                    "Master or equivalent (ISCED 7)",
                                    "Doctoral or equivalent (ISCED 8)"), 
          to = c("no-edu",
                 "Primary",
                 "Secondary",
                 "Secondary",
                 "Secondary",
                 "Tertiary",
                 "Tertiary",
                 "Tertiary",
                 "Tertiary"
                 ))

Create a new data frame for the model

data_sub <- subset.data.frame(data, select = c('D_INTERVIEW','dem3','dem5','dem6','dem8','dem9','H_URBRURAL', "edu"))

Create a new column as a dependent variable

data_sub['dem_val']<-(
  (as.numeric(data_sub$dem3)*sum(as.numeric(data_sub$dem3), na.rm = TRUE)  / 69867 /10) -
  (as.numeric(data_sub$dem5)*sum(as.numeric(data_sub$dem5), na.rm = TRUE)  / 69867 /10) + 
   (as.numeric(data_sub$dem6)*sum(as.numeric(data_sub$dem6), na.rm = TRUE)  / 69867 /10) - 
    
   (as.numeric(data_sub$dem8)*sum(as.numeric(data_sub$dem8), na.rm = TRUE)  / 69867 /10) + 
  (as.numeric(data_sub$dem9)*sum(as.numeric(data_sub$dem9), na.rm = TRUE)  / 69867 /10)) /5

Create a new data frame for the model

d_m <- data_sub[,c('edu',"dem_val", "H_URBRURAL")]

Delete missing values

d_m <- na.omit(d_m)

Histogram for dependent variable

hist(d_m$dem_val, main = "Histogram", ylab = "Anzahl", xlab = "democracy perception", col = "lightblue")

Education

d_m %>%
  group_by(edu) %>%
  ggplot(aes(dem_val, color=edu)) + 
  geom_histogram(aes(fill = edu), bins = 20) +
  facet_wrap(~edu) +
  theme_grey()+
  labs(x= "democracy perception",y = "count" )

Location

d_m %>%
  group_by(H_URBRURAL) %>%
  ggplot(aes(dem_val, color=H_URBRURAL)) + 
  geom_histogram(aes(fill = H_URBRURAL), bins = 20) +
  facet_wrap(~H_URBRURAL) +
  theme_grey()+
  labs(x= "democracy perception",y = "count" )

ggplot(d_m, aes(x=dem_val))+
geom_histogram(bins = 25)+
theme_grey()+
facet_wrap(edu~H_URBRURAL, ncol = 2)

Boxplot

ggplot(d_m, aes(H_URBRURAL, dem_val, fill=factor(edu))) +
geom_boxplot() +
theme_minimal()+
labs(fill = "Education", x="Location", y="Count")

d_m %>%
group_by(edu) %>%
  dplyr::summarize(Anzahl = n(), Mittelwert = mean(dem_val), Median = median(dem_val), Standardabweichung = sd(dem_val)) %>%
  mutate_if(is.numeric, round, 2)

At the non-education level, the average perception of democracy is 4.94 (SD = 1.47, n = 2967), At the primary level, the average perception of democracy is 5.03 (SD = 1.55, n = 8063), At the secondary level the perception is of democracy on average 5.39 (SD = 1.72, n = 31024), at the tertiary level of education the perception of democracy is on average 5.91 (SD = 1.78, n = 19122)

d_m %>%
group_by( H_URBRURAL) %>%
  dplyr::summarize(Anzahl = n(), Mittelwert = mean(dem_val), Median = median(dem_val), Standardabweichung = sd(dem_val)) %>%
  mutate_if(is.numeric, round, 2)

Beim Urban ist die Wahrnehmung der Demokratie im Durchschnitt 5.66 (SD = 1.81, n =40672), Beim Rural ist die Wahrnehmung der Demokratie im Durchschnitt 5.13 (SD = 1.53, n=20504).

d_m %>%
group_by(edu, H_URBRURAL) %>%
  dplyr::summarize(Anzahl = n(), Mittelwert = mean(dem_val), Median = median(dem_val), Standardabweichung = sd(dem_val)) %>%
  mutate_if(is.numeric, round, 2)

In urban the perception of democracy is on average 5.66 (SD = 1.81, n = 40672), in rural the perception of democracy is on average 5.13 (SD = 1.53, n = 20504).

Not; Although the difference seems small, it is tolerable since there is a large data set.

Profildiagramm

ggplot(d_m, aes(x=edu, y=dem_val, group=H_URBRURAL, color = H_URBRURAL))+
  stat_summary(fun.y = mean, geom="point", size=3)+
  stat_summary(fun.y = mean, geom="line")+
  stat_summary(fun.data = mean_cl_normal, geom="errorbar",width=.2, size=.25)+
  labs(x="Education", y="democracy perception")+
  theme_classic()

`fun.y` is deprecated. Use `fun` instead.`fun.y` is deprecated. Use `fun` instead.

Levene -Test

leveneTest(dem_val ~ edu*H_URBRURAL, data = data_sub, center = "mean")

In the present example, the Levene test is significant (F (7.61168) = 178.7, p = .000), so that heterogeneity of variance can be assumed. Since the variances are unfortunately not the same, it is recommended to carry out a correction using the Welch test.

without Welch-Korrektur

mehrAnova1 <- lm(dem_val ~ edu*H_URBRURAL, data = data_sub)
myAnova <- Anova(mehrAnova1 , type = 3)
myAnova

Anova Table (Type III tests)

Response: dem_val
               Sum Sq    Df    F value    Pr(>F)    
(Intercept)     29090     1 10106.2784 < 2.2e-16 ***
edu              3734     3   432.3802 < 2.2e-16 ***
H_URBRURAL         19     1     6.7261  0.009503 ** 
edu:H_URBRURAL    136     3    15.7329 3.181e-10 ***
Residuals      176066 61168                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Gesamtmodel wird signifikant.

with Welch-Korrektur

educ <- oneway.test(dem_val ~ edu, data = d_m, var.equal = F)
educ

    One-way analysis of means (not assuming equal variances)

data:  dem_val and edu
F = 756.52, num df = 3, denom df = 11969, p-value < 2.2e-16

There is a major effect that education levels have on the perception of democracy. (F (3.11969) = 756.52, p = .000). This means that the perception of democracy depends on the level of education

loc <- oneway.test(dem_val ~ H_URBRURAL, data = d_m, var.equal = F)
loc

    One-way analysis of means (not assuming equal variances)

data:  dem_val and H_URBRURAL
F = 1422, num df = 1, denom df = 47683, p-value < 2.2e-16

There is one major effect of the location on the perception of democracy. (F (1.47683) = 1422, p = .000). That means that the perception of democracy is dependent on the location

lage <- oneway.test(dem_val ~ edu*H_URBRURAL, data = d_m, var.equal = F)
lage

The interaction term of educational level and location on the perception of democracy is also significant (F (7,9917.6) = 436.4, p = .000).

Post-Hoc-Test

PostHoc <- aov(dem_val ~ edu*H_URBRURAL, data = d_m) 
TukeyHSD(PostHoc)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = dem_val ~ edu * H_URBRURAL, data = d_m)

$edu
                         diff          lwr       upr     p adj
Primary-no-edu     0.09201336 -0.001575966 0.1856027 0.0560253
Secondary-no-edu   0.45172974  0.367972963 0.5354865 0.0000000
Tertiary-no-edu    0.97206267  0.886060723 1.0580646 0.0000000
Secondary-Primary  0.35971638  0.305232930 0.4141998 0.0000000
Tertiary-Primary   0.88004931  0.822173726 0.9379249 0.0000000
Tertiary-Secondary 0.52033292  0.480260255 0.5604056 0.0000000

$H_URBRURAL
                  diff        lwr        upr p adj
Rural-Urban -0.3576343 -0.3861148 -0.3291538     0

$`edu:H_URBRURAL`
                                       diff         lwr         upr     p adj
Primary:Urban-no-edu:Urban       0.12013944 -0.05301552  0.29329440 0.4128865
Secondary:Urban-no-edu:Urban     0.47991700  0.32373712  0.63609689 0.0000000
Tertiary:Urban-no-edu:Urban      0.97526003  0.81772372  1.13279633 0.0000000
no-edu:Rural-no-edu:Urban       -0.16593574 -0.35985751  0.02798603 0.1581096
Primary:Rural-no-edu:Urban      -0.12800111 -0.29932414  0.04332193 0.3134280
Secondary:Rural-no-edu:Urban     0.10160783 -0.05830399  0.26151966 0.5329235
Tertiary:Rural-no-edu:Urban      0.44486775  0.27145570  0.61827979 0.0000000
Secondary:Urban-Primary:Urban    0.35977756  0.26928771  0.45026741 0.0000000
Tertiary:Urban-Primary:Urban     0.85512058  0.76230924  0.94793193 0.0000000
no-edu:Rural-Primary:Urban      -0.28607518 -0.43236910 -0.13978126 0.0000001
Primary:Rural-Primary:Urban     -0.24814055 -0.36280458 -0.13347652 0.0000000
Secondary:Rural-Primary:Urban   -0.01853161 -0.11532029  0.07825707 0.9991053
Tertiary:Rural-Primary:Urban     0.32472830  0.20696586  0.44249075 0.0000000
Tertiary:Urban-Secondary:Urban   0.49534302  0.44035702  0.55032902 0.0000000
no-edu:Rural-Secondary:Urban    -0.64585274 -0.77159623 -0.52010925 0.0000000
Primary:Rural-Secondary:Urban   -0.60791811 -0.69485116 -0.52098506 0.0000000
Secondary:Rural-Secondary:Urban -0.37830917 -0.43977106 -0.31684728 0.0000000
Tertiary:Rural-Secondary:Urban  -0.03504926 -0.12603007  0.05593155 0.9411005
no-edu:Rural-Tertiary:Urban     -1.14119577 -1.26862009 -1.01377144 0.0000000
Primary:Rural-Tertiary:Urban    -1.10326113 -1.19260814 -1.01391412 0.0000000
Secondary:Rural-Tertiary:Urban  -0.87365219 -0.93848349 -0.80882090 0.0000000
Tertiary:Rural-Tertiary:Urban   -0.53039228 -0.62368237 -0.43710219 0.0000000
Primary:Rural-no-edu:Rural       0.03793463 -0.10618633  0.18205560 0.9932806
Secondary:Rural-no-edu:Rural     0.26754357  0.13719381  0.39789334 0.0000000
Tertiary:Rural-no-edu:Rural      0.61080348  0.46420538  0.75740159 0.0000000
Secondary:Rural-Primary:Rural    0.22960894  0.13613706  0.32308082 0.0000000
Tertiary:Rural-Primary:Rural     0.57286885  0.45781697  0.68792074 0.0000000
Tertiary:Rural-Secondary:Rural   0.34325991  0.24601206  0.44050776 0.0000000

Since there is heterogeneity of variance, we should use Games Howell’s test. For this we need new columns as an interaction.

#first way
 f1 <- function(a,b){
  if (a == "Primary" & b == "Urban"){
   return("Primary:Urban")
  }else if(a == "no-edu" & b == "Urban"){
   return("no-edu:Urban")
  }else if(a == "Secondary" & b == "Urban"){
   return("Secondary:Urban")
  }else if(a == "Tertiary" & b == "Urban"){
   return("Tertiary:Urban")
  }else if (a == "Primary" & b == "Rural"){
   return("Primary:Rural")
  }else if(a == "no-edu" & b == "Rural"){
   return("no-edu:Rural")
  }else if(a == "Secondary" & b == "Rural"){
   return("Secondary:Rural")
  }else if(a == "Tertiary" & b == "Rural"){
   return("Tertiary:Rural")
  }
 }

d_m['eu'] <- mapply(f1, d_m$edu, d_m$H_URBRURAL)

#second way
d_m$eu <- paste(d_m$edu,d_m$H_URBRURAL)

head(d_m)

gm <- rbind(games_howell_test(d_m, dem_val ~ edu),games_howell_test(d_m, dem_val ~ H_URBRURAL),games_howell_test(d_m, dem_val ~ eu))
gm

subset(gm, subset = p.adj.signif == "ns")

subset(gm, subset = estimate == max(gm$estimate))

eta <- effectsize::eta_squared(mehrAnova1, partial = TRUE)
eta

Parameter      | Eta2 (partial) |       90% CI
----------------------------------------------
edu            |           0.03 | [0.03, 0.04]
H_URBRURAL     |           0.01 | [0.01, 0.01]
edu:H_URBRURAL |       7.71e-04 | [0.00, 0.00]

for (i in 1:length(eta$Parameter)){
  st <- sqrt(eta$Eta2_partial[i]/(1-eta$Eta2_partial[i]))

  a <- sprintf("Effect size f=%s fürs f= %.3f", eta$Parameter[i],st)
  print(a)
}

[1] "Effect size f=edu fürs f= 0.189"
[1] "Effect size f=H_URBRURAL fürs f= 0.103"
[1] "Effect size f=edu:H_URBRURAL fürs f= 0.028"

It turns out that there is a difference between the level of education and the perception of democracy (F (3.11969) = 756.52, p <.000). H1 for the main effect A is accepted. In addition, depending on the location, there is a difference between the perception of democracy (F (1.47683) = 1422, p <.000). H1 for the main effect B is accepted. The interaction term of educational level and situation on the perception of democracy is significant (F (7,9917.6) = 436.4, p <.000). Therefore H1 for the interaction effect AXB is accepted

The PostHoc was carried out with Games Howell. It shows that all groups (with the exception of Primary and non-education (p = 2.10e-02)) differ significantly in terms of educational level and location. And there was a significant difference between primary and non-education. In relation to the Games Howell table of results, the differences between the combination of non-education and primary with location also appear to be insignificant. In general,the differences between the Urban and the Rural is also based on education level, and rural areas follows urban areas one step behind. For instance, the differences between Urban with secondary level education and Rural with tertiary level is not clear.Finally, as expected, the biggest difference is between non-education rural and tertiary urban (e = 1.14, p = 0)

The effect sizes are a weak effect for the main effect of educational level and location (f = 0.189, f = 0.103)